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Recovery  of  Eurasian  Crustal  and  Upper-Mantle  Structure  by 
Higher-Mode  Waveform  Analysis 

Thomas  H.  Jordan,  Arthur  l.  Lerner-Lam*.  and  Lind  S.  gee 

Department  of  Earth,  Atmospheric  and  Planetary  Sciences 
Massachusetts  Institute  of  Technology,  Cambridge,  MA  02139 


This  paper  reviews  the  waveform  inversion  technique  of  I^rner-Lam  and  Jordan  [1983]  and  reports  on  its 
application  to  the  study  of  Eurasian  crustal  and  upper-mantle  structure.  We  have  successfully  fit  complex 
P-SV  wavegroups  on  vertical-component  seismograms  for  paths  crossing  northern  Eurasia  with  a  model 
having  a  39-km  crust  and  no  asthenospheric  low-velocity  zone.  Between  the  Moho  and  the  400-km 
discontinuity,  the  shear  velocities  found  from  the  P-SV  waveform  analysis  are  consistently  lower  than 
those  inferred  from  the  SH  waveform  modeling  of  Grand  and  Helmberger.  We  suggest  that  this  discrepancy 
is  diagnostic  of  a  polarization  anisotropy  associated  with  the  olivine-rich  mineralogy  of  the  thick, 
basalt-depleted  chemical  boundary  layer  that  characterizes  the  upper  mantle  beneath  stable  continents,  a 
hypothesis  consistent  with  the  higher-mode  results  of  Cara  and  others.  The  algorithms  employed  in  this 
study  can  be  generalized  to  incorporate  such  anisotropic  structures  and  can  be  applied  to  the  construction  of 
fully  three-dimensional  models  of  Eurasia,  as  well  as  to  the  study  of  upper-mantle  attenuation  structure. 


INTRODUCTION 

Seismological  methods  for  the  detection  and  discrimination  of  underground  nuclear  explosions 
rely  on  models  of  the  crustal  and  upper-mantle  structure,  particularly  attenuation  structure  in  the 
vicinity  of  the  testing  areas.  Estimates  of  the  body-wave  magnitude  bias  between  the  Eastern  Kazakh 
Test  Site,  where  the  recent  large  Soviet  shots  have  been  fired,  and  the  Nevada  Test  Site,  which 
provides  the  bulk  of  the  magnitude-yield  calibration  data,  range  from  0.2  to  0.4  magnitude  units, 
corresponding  to  factors  of  about  1.6  to  2.8  in  yield.  The  uncertainty  in  this  bias  dominates  the 
errors  in  the  yields  determined  from  teleseismic  P  waves  and  has  proven  to  be  a  source  of 
controversy  in  verifying  compliance  with  the  Threshold  Test  Ban  Treaty. 

For  fixed  receiver  networks,  the  body-wave  attenuation  bias  is  primarily  controlled  by  local 
and  regional  variations  in  upper-mantle  structure  beneath  the  sources.  These  variations  appear  to  be 
largest  in  the  asthenosphere  below  100  km.  Recent  models  of  the  subcontinental  mantle  suggest  that 
substantial  variations  in  temperature  and  composition  extend  to  depths  exceeding  200  km  and  that 
these  variations  correlate  with  the  long-term  tectonic  history  of  the  overlying  crust  [Jordan,  1981]. 
Therefore,  the  problems  of  yield  estimation  are  coupled  to  fundamental  questions  regarding 
continental  evolution. 

Previous  work  on  mapping  the  three-dimensional  structure  of  Eurasia  has  employed 
fundamental-mode  surface  v/aves  [e.g.,  Feng  and  Teng,  1983],  which  lack  sufficient  resolution 
below  100-km  depth,  or  a  limited  number  of  multiply  reflected  body  waves  with  restricted 
geographical  coverage  [e.g.,  Sipkin  and  Jordan,  1976;  Burdick  et  ai,  1983].  Grand  et  al.  [1985] 
have  recently  extended  the  SS-S  technique  of  Burdick  et  al.  [1983]  and  Grand  and  Helmberger 
[1984]  to  include  SSS  phases  and  thereby  increase  the  geographical  range  of  the  body-wave 
observations,  but  their  work  has  thus  far  been  confined  to  the  forward  modeling  of  SH-polarized 
signals. 


*  Now  at  Lament- Doherty  Geological  Observatory,  Palisades,  NY  10964 


*,  A 


.•vy-.-' 


In  a  DARPA/AFGL-sponsored  project  initiated  this  year,  we  are  attempting  to  develop 
techniques  for  the  systematic  inversion  of  both  P-SV  and  SH  waveforms  from  higher-mode 
dispersed  wave  groups  (which  include  multiply  reflected  body  waves)  for  lateral  heterogeneity.  These 
methods  are  capable  of  giving  much  better  vertical  resolution  of  path-averaged  properties  than  the 
dispersion  of  the  fundamental-mode  surface  waves  alone.  They  will  be  applied  to  large  data  sets 
collected  from  GDSN,  NARS,  WWSSN  and  other  seismic  networks  to  obtain  three-dimensional 
models  of  Eurasian  crustal  and  upper-mantle  structure.  Here  we  describe  the  formulation  of  the 
waveform  inversion  procedure  and  give  the  results  of  some  preliminary  modeling  of  the  structure 
along  paths  traversing  the  northern  part  of  Eurasia. 


FORMULATION  OF  THE  WAVEFORM-INVERSION  TECHNIQUE 

Our  formulation  is  based  on  the  mode-isolation  and  waveform-inversion  techniques  described 
by  Lerner-Lam  and  Jordan  [1983].  The  observed  seismogram  is  represented  as  the  sum  of 
fundamental  (n  =  0)  and  higher-mode  (n  ^  1)  surface  waves: 

=  X 

n=Q 

u  Jt)  is  the  seismogram  for  the  n*  mode  branch.  A  synthetic  seismogram  s  (t)  is  calculated  for  a 
cnosen,  spherically  symmetric  earth  model  m  (r),  and  the  difference  between  the  observed  and  the 
synthetic  is  formed: 

N 

As(t)  s  s(t)  -  i(r)  =  Am„(0 

n=0 

A«^(r)  is  the  differential  seismogram  for  the  n*  mode  branch. 

To  increase  the  signal-to-noise  ratio  (SNR)  for  a  specified  mode  and  reduce  the  interference 
from  spurious  signals  and  other  modes,  we  construct  matched  filters  from  the  synthetic  branch 
seismograms.  At  the  lag  time  x  we  define  the  observed  branch  cross-correlation  function  (BCCF) 


-oo 

UmiO  s{t  + dt 

^  -OO 


and  the  synthetic  BCCF 


provides  an  approximate  description  of  the  mode-mode  interference,  so  an  appropriate  data  functional 
for  the  structural  inverse  problem  is  the  differential  BCCF, 


In  our  procedure,  we  seek  a  model  perturbation  which  to  first  order  minimizes  a  quadratic  form  in 
The  quadratic  form  includes  a  symmetric  taper  about  zero  lag,  so  the  points  near  r  =  0,  where 
the  SNR  is  greatest,  receive  the  most  weight 

For  the  first-orbit,  vertical-component  surface  waves  received  at  stations  far  away  from  the 
source  and  its  antipode,  branch  seismograms  can  be  approximated  by  an  integral  over  continuous 
wavenumber  A,  whose  asymptotic  form  is 


(A,0  =  G„  (A,A,0  cos  (A,A,/)  dk 
•'0 

A  is  epicentral  distance,  and  and  are  the  amplitude  and  phase  kernels 


G„(X,A,0  =  - ;  exp[-a„(X)f] 

27C  (sinA)*^ 


^^(A,A,0  =  AA  -  71/4  -  ©„(A)r  +  <I>„ 


and  are  displacement  and  excitation  scalars,  is  the  source  phase,  is  the  dispersion 
function,  and  is  the  decay  function.  and  depend  on  the  source  mechanism, 

whereas  and  depend  on  the  path-averaged  elastic  and  anelastic  parameters,  respectively.  The 
horizontal-component  seismograms  have  similar  forms. 

We  have  developed  an  efficient  algorithm  for  computing  branch  synthetics  which  is  based  on 
Filon  quadrature.  The  algorithm  is  adaptive  and  reduces  computation  time  by  an  order  of  magnitude 
over  conventional  mode-summation  techniques  \Lerner-Lam  andJordan,  1983]. 

Departures  of  the  real  earth  from  the  spherically  symmetric  reference  model  m  are  represented 
as  perturbations  to  the  amplitude  and  phase  kernels: 

G„(A,A,r)  =  G„(A,A,0  [1+Y„(A)] 

=  4„(X,A,f)  -  A©„(A)t 

Here  y  (A)  is  the  perturbation  to  the  relative  amplitude  of  the  mode,  and  A©^  =  ©  -  w„  is  the 
perturoation  to  its  dispersion  function.  If  the  perturbations  to  the  displacement  scalar  and  the 
excitation  scalar  can  be  ignored,  a  good  approximation  when  the  lateral  variations  along  the  path 
are  small,  then  the  relative  amplimde  perturbation  can  be  written 

Yn(^)  =  -Aa„(A)r 

where  Aa„  =  %  -  Aa^,  which  in  turn  can  be  related  to  the  perturbation  in  the  specific  attenuation 
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(X,).  In  terms  of  the  differential  BCCF,  the  linearized  forward  problem  becomes 


I 


n=0 


c°° 

•^0 


The  partial  derivatives  and  can  be  evaluated  by  the  same  adaptive  Filon  algorithm  used  to 
compute  the  synthetic  seismograms.  The  integrals  over  wave-number  are  thereby  discretized  on  the 
Filon  grids  {X^:  1  =  1,  2, . . .  ,L^}.  To  set  up  a  matrix  system  describing  the  differential 
BCCFs  for  M  +  l  mode  branches  (m  =  0,  1, .  .  . ,  M  ^  N)  at  each  of  P  stations  (jj  =  1, 
1, ...  ,P),  we  discretize  the  cross-correlations  at  the  set  of  lag  times  {t^  =  it  At:  k  = 
-K,  -AT+l, .  .  .  ,  0  ,  .  .  . ,  K-l,  K}  and  define  the  vectors 

[AS]^,  =  AS^(x,) 

[Acd]„,  =  A0)„(X.,) 


and  the  matrices 

^^^mpk,nJl  ~  ^mnp  ^1^ 

^^^ntpk,nA  ~  ^mnp  '^0 

This  yields  the  linearized,  discretized  forward  problem 


AS  =  BA©  +  DAa 


(1) 


The  matrices  appearing  in  this  linear  system  have  I  =  (M+\)-P  -{2K+\)  rows  and 
N 

J 

n=0 


columns.  For  typical  ranges  of  indices  used  in  our  work  to  date,  these  dimensions  are  on  the  order 
of  10"*  and  10^  ,  respectively. 

The  vector  AS  can  be  computed  from  a  set  of  observed  seismograms,  and  equation  (1)  inverted 
for  the  perturbations  to  the  dispersion  and  the  attenuation  functions.  In  practice,  it  is  convenient  to 
parameterize  the  vectors  A©  and  Aa  by  perturbations  Am  to  the  radial  starting  model.  This  reduces 
the  dimension  of  the  system  and  constrains  the  perturbed  seismograms  to  correspond  to  a  realistic 
earth  structure.  We  interpret  the  resulting  one-dimensional  model  as  the  average  of  the 
three-dimensional  structure  ^ong  the  path  between  the  source  and  receiver.  Inversion  procedures 
based  on  this  representation  have  been  dscussed  by  Lerner-Lam  and  Jordan  [1983]. 
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Application  To  Eurasian  Data 

The  mode-isolation  and  waveform-inversion  technique  has  been  applied  to  seismograms 
recorded  at  five  WWSSN  stations  in  western  Europe  from  four  earthquakes  in  the  Kuril-Kamchatka 
Seismic  Zone  (Figure  1).  The  focal  depths  for  the  events  range  from  134  km  to  544  km,  and  the 
epicentral  distances  from  55°  to  78°.  The  data  set  comprises  a  total  of  fourteen  paths  traversing  a 
variety  of  tectonic  structures  in  northern  Eurasia,  from  the  tectonically  active  region  east  of  the 
Verkhoyansk  Suture  Zone  to  the  Baltic  Shield.  The  long-period,  vertical-component  seismograms 
were  digitized,  response-normalized,  and  low-pass  filtered  with  a  comer  at  35  mHz. 


Figure  1.  Azimuthal  equidistant  projection  of  the  event-station  disuibution  for 
the  northern  Eurasian  path,  with  the  pole  centered  on  the  epicenter  of  Event  1.  The 
events,  shown  as  triangles,  range  in  focal  depth  from  134  to  544  km.  The  stations 
are  shown  as  circles. 


In  our  original  analysis  of  this  data  set  [Lerner-Lam  and  Jordan,  1983],  we  chose  the  very  smooth 
structure  of  Cara  et  al.  [1980]  as  a  starting  model  and  solved  for  shear- velocity  perturbations  that 
were  a  smooth  function  of  depth.  The  modeling  presented  in  this  report  is  based  on  a  layered 
structure  with  four  upper-mantle  discontinuities  below  the  crust-mantle  boundary  and  polynomial 
variations  in  the  seismic  velocities  and  density  of  each  layer;  the  polynomials  were  taken  to  be  linear 
throughout  the  upper  mantle.  Perturbations  were  allowed  in  the  depths  of  all  discontinuities,  as  well 
as  in  the  parameters  of  the  polynomials.  The  starting  model  was  set  equal  to  PREM  [Dziewonski  and 
Anderson,  1983]  in  the  lower  mantle  and  conformed  to  regionalized  body-wave  and  surface-wave 
models  of  continental  structure  above  this  depth;  it  contained  no  low  velocity  zone. 

The  attenuation  model  used  in  our  calculations  was  that  of  Masters  and  Gilbert  [1983].  In 
these  preliminary  experiments,  we  did  not  attempt  to  invert  for  attenuation  structure.  We  did, 
however,  desensitize  ^e  inversion  to  possible  amplitude  fluctuations  associated  with  Q  variations  and 
source  effects  using  the  normalization  scheme  described  by  Lerner-Lam  and  Jordan  [1983,  sect.  5.2], 
which  employs  a  projection  operator  to  annihilate  the  second  term  on  the  right-hand  side  of  equation 
(1)  for  a  restricted  class  of  amplitude  perturbations. 

The  model  obtained  from  waveform  inversion,  EU2,  is  plotted  in  Figure  2.  Comparisons  of 
observed  and  synthetic  waveforms  and  BCCFs  are  displayed  in  Figures  3-6.  The  structure  is  a 
simple  one,  but  it  provides  a  remarkably  good  fit  to  the  phase  of  the  waveforms  over  the  entire  range 
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of  focal  depths,  including  both  the  fundamental  Rayleigh  mode  and  the  complex  wave  groups 
respresenting  the  interference  of  the  higher  modes.  The  pulses  seen  on  the  seismograms  in  Figure  6, 
for  example,  where  the  match  is  nearly  wiggle-for-wiggle,  are  dominated  by  contributions  from  the 
third  and  fourth  overtones.  These  waveforms  can  also  be  represented  as  the  sum  of  multiply 
reflected  P-SV  body  phases,  predominantly  SS  and  SSS;  the  agreement  shows  that  the  model 
satisfies  their  travel  times.  Our  waveform  analysis  is  thus  the  P-SV  equivalent  of  the  SH  analysis 
performed  by  Grand  and  Helmberger  [1984a,b] ,  and  our  resolution  of  sub-lithospheric  structure 
should  be  comparable  to  theirs.  (Because  our  technique  is  formulated  as  a  linearized  inverse 
problem,  its  resolving  power  can  be  analyzed  using  form^  linear  methods,  which  is  one  of  its  major 
advantages  over  forward-modeling.  Preliminary  cdculations  have  been  made  by  Lerner-Lam  [1983], 
and  a  more  detailed  resolving-power  study  is  given  by  Gee  et  al..  [1985].) 


EU2 

VELOCITY  (km  s'’).  DENSITY  (gm  cm'*) 


Figure  2.  EU2,  a  path-averaged  model  of  northern  Eurasia  obtained  by  waveform 
inversion  of  fundamental  and  higher-mode  surface  waves. 


EU2  shear  velocities  are  compared  with  the  SNA  model  of  Grand  and  Helmberger  [1984a]  in  Figure 
7.  The  latter  structure  was  derived  from  S  and  SS  waves  traversing  the  Canadian  Shield,  but  Grand 
etal.  [1985]  have  shown  that  it  also  provides  a  good  description  of  SH  propagation  across  the  stable 
platforms  of  Eurasia.  The  two  models  have  comparable  velocity  jumps  at  400  km  and  are  in  good 
agreement  below  this  depth.  However,  above  400  km,  there  are  significant  differences  between  the 
two  structures.  SNA  has  a  low-velocity  zone  centered  at  about  200-km  depth,  whereas  the  shear 
velocities  in  EU2  increase  monotonically  throughout  the  upper  mantle.  Experiments  with  this  and 
other  parameterizations,  including  structures  parameterized  by  continuous  variations  with  depth, 
indicate  that  our  data  set  cannot  be  satisfied  by  a  model  with  an  LVZ  as  pronounced  as  SNA;  the 
velocity  variation  between  100  km  and  200  km  is  tightly  constrained  by  the  difference  in  the  arrival 
times  of  the  fundamental  and  first  higher-mode  wave  groups,  which  are  well  observed  on  many  of 
the  seismograms  used  in  this  study. 

The  most  obvious  discrepancy  between  the  two  models  is  the  negative  offset  in  the  shear 
velocity  of  EU2  relative  to  SNA  throughout  the  mantle  above  400  km,  which  averages  0.16  km/s. 
The  differences  are  largest  in  the  uppermost  mantle  (the  lid  of  SNA),  where  they  locally  reach  0.27 
km/s,  but  remain  substantial  even  below  200  km.  Given  the  high  resolution  of  both  methods  in  the 
interval  200-400  km,  the  differences  appear  to  be  real.  For  example,  inversions  of  our  data  set  with 
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Figure  3.  Comparisons  of  observed  and 
syn-hetic  waveforms  (top  panel)  and  BCCFs 
(bottom  panel)  recorded  at  NUR  for  Event  1 
(h  =  134  km).  BCCFs  plotted  a.s  solid  lines 
observed,  those  plotted  as  dotted  traces  are 
synthetic.  The  BCCF  mode  number  is  to  the 
right  of  each  pair  of  traces  in  the  lower  panel 
The  energy  in  the  waveforms  is  concentrated  in 
the  fundamental  mode,  which  appears  as  the 
well-dispersed  wavetrain  beginning  at  about 
29  min,  although  the  S  and  SS  arrivals  are  also 
visible  at  20  min  and  25  min. 
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Figure  4.  Comparisons  of  observed  and 
synthetic  waveforms  and  BCCFs  for  Event  2 
(h  =  181  km)  recorded  at  KON.  Increased 
relative  amplitude  levels  of  the  higher-modes 
are  due  to  the  increase  in  source  depth. 
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Figure  5.  Comparisons  of  observed  and 
synthetic  waveforms  and  BCCFs  for  Event  3 
(h  =  344  km)  recorded  at  UME.  Fundamental 
mode  energy  is  still  visible  at  29  minutes,  but 
the  record  is  dominated  by  higher-mode  arrivals. 
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Figure  6.  Comparisons  of  observed  and 
synthetic  waveforms  and  BCCFs  for  Event  4 
(h  =  544  km)  recorded  at  NUR.  Fundamental 
mode  enery  is  almost  totally  absent  from  the 
waveforms;  complex  higher-mode  arrivals 
dominate  the  signal.  Near  wiggle-for- wiggle 
fits  are  seen  for  the  second,  third,  and  fourth 
higher  modes. 
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the  velocities  in  this  layer  constrained  at  SNA  values  yield  very  poor  fits  to  the  observed 
seismograms. 

&me  of  the  discrepancy  may  result  from  the  lower  velocities  encountered  along  the  portions  of 
the  surface-wave  paths  fraversing  the  tectonically  active  areas  of  northeastern  Asia,  but  we  are 
skeptical  that  an  explanation  based  on  path  differences  can  account  for  it  entirely.  Grand  ct  al.  [  1985J 
have  modeled  the  transition  from  the  active  foldbelts  of  central  Asia  to  the  Russian  platform  using 
their  multiple-S  technique,  and  the  loM/est  average  velocities  they  find  are  still  higher  than  those  of 
EU2,  even  though  the  latter  is  derived  from  a  data  set  primarily  sampling  the  high-velocity  shields 
and  platforms  of  northern  Eurasia. 
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Figure  7.  Comparison  of  EU2  shear  velocities  (solid  line)  with  model  SNA 
(dotted  line)  of  Grand  and  Helrrtberger  (1984a).  SNA  was  derived  from  SH -polarized 
S  and  SS  waves  traversing  the  Canadian  Shield,  but  is  representative  of  SH 
propagation  across  stable  Eurasian  platforms  {Grand  el  al.  1985).  EU2,  derived  from 
vertical-component  higher-mode  surface  waves,  is  representative  of  P-SV 
propagation.  The  structures  are  virtually  identical  below  the  400  km  discontinuity, 
but  exhibit  distinct  velocity  differences  in  the  upper  mantle  above  400  km. 

Although  a  definitive  statement  must  await  the  treatment  of  both  P-SV  and  SH  data  sets  from 
common  paths  by  the  same  inversion  technique  (which  is  one  of  our  goals  for  the  next  year),  we 
suspect  that  the  discrepancy  in  the  average  shear  velocity  between  EU2  and  SNA  is  diagnostic  of 
deep-seated  polarization  anisotropy  in  the  Eurasian  upper  mantle.  This  hypothesis  is  consistent  with 
the  observations  and  modeling  results  of  Cara  etal.  [1980],  who  found  similar  differences  in  SV  and 
SH  velocities  from  the  inversion  of  higher-mode  Rayleigh  and  Love  waves.  We  speculate  that  this 
anisotropy  is  associated  with  the  existence  of  a  thick,  basalt-depleted  (and  therefore  olivine-rich) 
chemical  boundary  layer  beneath  Eurasia  [Jordan,  1981], 
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Future  work 

The  waveform-inversion  technique  described  in  this  report  will  be  extended  to  include 
polarization  anisotropy  and  variations  in  attenuation  structure.  We  will  apply  the  method  to 
three-component  data  collected  from  GDSN,  NARS,  WWSSN  and  other  seismic  networks  to  obtain 
models  of  the  crust  and  upper  mantle  over  various  Eurasian  paths.  These  models  will  be  employed  to 
constrain  the  three-dimensional  structure  of  the  Eurasian  continent 
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In  three-dimensional  studies  of  earth  structure  based  on  surface-wave  tomography,  horizontal 
resolution  is  limited  primarily  by  the  distribution  of  paths,  while  vertical  resolution  is  governed  by  the 
distribution  of  mode  types.  To  obtain  good  vertical  resolution  below  200-km  depth  requires  the  use  of 
higher  modes.  This  paper  investigates  the  vertical  resolving  power  of  source  and  receiver  arrays  in 
determining  Eurasian  crustal  and  upper-mantle  structure  using  the  waveform  inversion  technique  of 
Lerner-Lam  and  Jordan  [1983].  Unlike  approaches  to  higher-mode  analysis  which  rely  on  stacking  to 
separate  modes  in  the  Irequency-wavenumber  domain,  this  procedure  isolates  the  modes  in  the  time  domain 
by  cross  correlating  a  mode-branch  synthetic  with  the  observed  and  synthetic  seismograms.  The  difference 
between  the  observed  and  synthetic  branch  cross-correlation  functions  is  approximated  as  a  linear  functional 
of  the  residual  dispersion,  which  is  parameterized  in  turn  by  perturbations  to  earth  structure.  These 
differential  branch  cross-correlation  functions  are  inverted  for  model  perturbations,  from  which  new 
dispersion  curves  and  synthetic  seismograms  are  computed,  and  the  process  is  iterated  until  convergence. 
This  formulation  permits  evaluation  of  real  and  hypothetical  data  sets  by  the  resolution  and  covariance 
analysis  of  solutions  to  the  linearized  inverse  problem.  The  results  show  that  good  vertical  resolution  of 
upper-mantle  structure  can  be  obtained  with  this  method  from  sparse  arrays  of  sources  and/or  receivers. 
Source  arrays  are  particularly  effective  in  enhancing  resolution,  provided  that  the  source  depths  are  well 
distributed  and  the  source  centroids  and  moment  tensors  are  well  determined.  Resolution  at  depth  may  be 
obtained  with  this  technique  even  with  shallow  sources,  as  long  as  the  higher-mode  branch  cross-correlation 
functions  are  weighted  appropriately. 


INTRODUCTION 

Studies  of  multiply-reflected  shear  waves,  such  as  ScS  [Sipkin  and  Jordan,  1976]  and  SS 
[Grand  et  ai,  1985],  as  well  as  studies  of  Rayleigh-wave  dispersion  [Feng  and  Teng,  1983], 
indicate  that  there  are  substantial  lateral  variations  in  upper-mantle  structure  underlying  the  Eurasian 
continent  to  depths  exceeding  200  km.  Recent  models  of  the  composition  and  development  of  the 
subcontinental  mantle  [Jordan,  1978,  1981;  Davies,  1979',  Richter,  1985]  attribute  this 
heterogeneity  to  compositional  and  thermal  differences  associated  with  the  structure  of  what  Jordan 
[1975]  has  termed  the  'continental  tectosphere.'  These  models  predict  a  strong  correlation  between 
surface  geological  features  and  upper-mantle  heterogeneity,  a  hypothesis  which  may  be  tested  by 
detailed  seismic  imaging  of  the  Eurasian  crust  and  upper  mantle. 

One  approach  to  such  three-dimensional  studies  in  Eurasia  employs  fundamental-mode 
Rayleigh  waves  to  constrain  lateral  variations  in  structure  [Feng  and  Teng,  1983].  While  good 
horizontal  resolution  is  possible  with  this  technique,  fundamental-mode  dispersion  is  not  sensitive 
to  details  in  vertical  structure  at  depths  greater  than  200  km.  In  principle,  the  addition  of  higher¬ 
mode  data  can  improve  structural  resolution,  but  mode-mode  interference  complicates  the 
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measurement  of  dispersion  in  both  the  lime  and  frequency  domain  [Cara,  1978;  Lerner-Lam. 
1982],  One  method  of  analysis  of  higher- mode  information  requires  sophisticated  processing  of 
the  data  in  the  frequency-wavenumber  domain  to  separate  the  modes  [Nolet,  1975,  1976;  Cara, 
1978,  1979;  Chou  and  Dziewonski,  1980].  While  this  technique  has  been  used  successfully  in 
studies  of  one-dimensional  structure  in  Eurasia  [Nolet,  1975;  Cara  et  al.,  1980;  Chou  and 
Dziewonski,  1981],  the  Pacific  [Cara,  1979],  and  North  America  [Cara,  1978,  1979],  the 
procedure  requires  an  array  of  stations  distributed  over  several  thousand  kilometers  and 
more-or-less  aligned  along  the  path  from  the  epicenter.  This  requirement  restricts  the  number  of 
paths  which  may  be  studied  with  existing  seismic  networks  and  therefore  limits  the  applicability  of 
the  method  to  tomographic  experiments. 

A  different  approach  to  higher-mode  analysis  has  been  taken  by  Lerner-Lam  and  Jordan 
[1983].  They  employ  waveform  inversion  to  model  the  difference  between  an  observed  and 
synthetic  seismogram  as  a  measure  of  the  departure  of  an  earth  model  from  an  average  of  structure 
along  the  path.  These  differential  seismograms  are  parameterized  by  perturbations  to  the 
fundamentd  and  higher-mode  dispersion,  which  are  relat^  to  perturbations  to  the  earth  model.  In 
traditional  waveform  inversion  methods,  the  differential  seismograms  are  inverted  directly  for 
model  perturbations,  with  some  weighting  scheme  applied  on  a  time  point  by  time  point  basis. 
However,  Lerner-Lam  and  Jordan  [1983]  isolate  higher-mode  information  explicitly  by  cross 
correlating  the  differential  seismograms  with  the  mode-branch  synthetics  and  windowing  the  result 
in  the  lag-time  domain.  In  addition  to  enhancing  the  signal-to-noise  ratio  for  a  particular  mode  and 
reducing  interference  from  spurious  signals  and  other  modes,  this  approach  permits  the 
contribution  from  each  mode  to  be  assessed  and  weighted  individually.  The  forward  problem 
relating  the  branch  cross-correlation  functions  to  model  perturbations  is  linearized  by  the  application 
of  first-order  perturbation  theory.  Solutions  are  derived  by  a  generalized  least-squares  inverse,  and 
the  process  is  iterated  to  convergence.  The  information  provided  by  this  approach  is  equivalent  to 
the  waveform  analysis  of  multiply-reflected  body  waves  [Grand  and  Helrnberger,  I984a,b;  Grand 
et  al.,  1985],  which  uses  trial-and-error  forward  modelling  techniques.  However,  the  ability  of 
these  methods  to  resolve  upper-mantle  structure  is  difficult  to  assess,  whereas  the  linearized 
formulation  used  in  this  study  permits  the  evaluation  of  the  solutions  by  the  now  standard  resolving 
power  analysis  of  linear  inverse  theory. 

In  the  next  few  years,  millions  of  dollars  will  be  spent  on  the  deployment  of  new  digital 
seismic  instruments  around  the  world.  In  particular,  the  Incorporated  Research  Institutes  for 
Seismology  (IRIS)  network  and  the  Portable  Array  of  Seismometers  for  Studies  of  the  Crust  and 
Lithosphere  (PASSCAL)  will  be  coming  on-line  in  the  near  future.  In  order  to  maximize  the 
benefit  from  this  expansion,  it  is  important  to  understand  how  various  configurations  of  sources 
and  receivers  will  affect  the  resolution  of  earth  structure.  For  example,  an  array  of  seismometers 
recently  has  been  established  in  Western  Europe.  The  Network  of  Autonomously  Registering 
Stations  (NARS)  array  spans  2500  km  with  14  receivers  [Nolet,  Dost,  and  Paulssen,  1984].  One 
question  which  needs  to  be  addressed  is  whether  dense,  linear  arrays  such  as  NARS  are  essential 
for  detailed  studies  of  upper-mantle  structure.  While  the  higher-mode  wavenumber-stacking 
techniques  of  Nolet  [1975,  1976]  and  Cara  [1978,  1979]  require  this  type  of  array  to  obtain  the 
necessary  spatial  coverage,  the  waveform  inversion  technique  of  Lerner-Lam  and  Jordan  [1983], 
which  makes  use  of  a  priori  information  about  structure  along  the  path,  does  not.  Consequently,  a 
more  efficient  deployment  of  such  receiver  arrays  would  distribute  them  over  a  larger  area  in  order 
to  sample  the  greatest  number  of  source-receiver  paths. 


FORMULATION  OF  THE  INVERSE  PROBLEM 

In  the  waveform  inversion  technique  of  Lerner-Lam  and  Jordan  [1983],  an  observed 
seismogram  is  characterized  as  the  sum  of  the  fundamental  {n  =  0)  and  higher-mode  {n  >  0) 


surface  waves: 


s(.T,t)  =  2)M„(r,0  (1) 

where  M^(r,r)  is  the  seismegram  for  the  n*  mode  branch  at  a  receiver  position  r  and  time  t.  A 
synthetic  seismogram,  ^  (r,f)  ,  is  calculated  for  a  particular  spherically-symmetric  earth 
model,  m(r),  complete  to  radial  oider  N: 

N  _ 

s(r,0  =  2  (2) 

n=0 

where  a^(r,f)  is  the  synthetic  seismogram  for  the  n*  mode  branch.  In  our  application  of  this 
technique,  the  sum  over  N  is  truncated  at  the  seventh  higher  mode  (Figure  1).  Although  phases  of 
high  apparent  velocity,  such  as  core  reflections,  are  not  well  represented  as  a  consequence,  these 
phases  are  nearly  transversely  polarized  and  do  not  contribute  to  the  vertical-component 
seismograms  used  in  this  study.  u„(r,0  is  formulated  using  an  asymptotic  approximation  to 
Gilbert's  [1976]  exact  travelling-wave  representation  of  the  displacement  field: 

S 

Mrt(r,t)  ~  5^  G^(^,r,r)  cos  ^  •  (3) 

Jq 

In  this  formulation,  X  is  the  surface  spherical  wavenumber,  and  G^(K,r.t)  and  ^(A.,r,r)  are  the 
amplitude  and  phase  kernels  for  the  orbit  of  the  higher  mode.  Since  we  are  concerned  with 
arrivals  corresponding  to  the  first  (minor-arc)  orbit,  we  t^e  s  =  1  and  drop  the  orbital  indices. 

The  difference  between  an  observed  and  synthetic  seismogram,  As(r,r),  can  be  approximated 
as  a  radial  perturbation,  Am(r),  to  the  starting  earth  model,  m(r),  where 

N 

As(r,f)  =  5(r,0  -  .<r,r)  =  ^Au^(r,t)  (4) 

n=0 

and  Au^(r,t)  is  the  differential  seismogram  for  the  n*  mode  branch.  In  order  to  estimate  Am(r) 
from  A5(r,t),  it  is  necessary  to  minimize  both  observational  and  representational  errors.  These 
include  random  seismic  noise  and  digitizing  errors,  as  well  as  error  processes  which  scale  with 
signal  strength,  such  as  inaccuracies  in  representing  the  source,  the  instrument  response,  and  wave 
propagation  through  the  earth.  One  method  of  enhancing  the  signal  corresponding  to  a  particular 
mode  and  reducing  the  effect  of  ambient  noise  and  spurious  signals  is  to  construct  matched  filters 
[Chou  and  Dziewonski,  1980]  by  cross  correlating  the  differential  seismograms  with  the 
mode-branch  synthetics  and  windowing  the  results  in  the  lag- time  domain.  Lerner-Lam  and  Jordan 
[1983]  define  ^e  observed  branch  cross-correlation  function. 


S^{T,x)  =  ujx,t)  *  ^(r,0 


(5) 


1  3 


and  the  synthetic  branch  cross-correlation  function, 


5^(r,x)  =  ujr,t)  *  sir.t)  (6) 

where  x  is  the  lag  time.  With  this  notation,  an  appropriate  functional  for  the  structural  inverse 
problem  is  the  differential  branch  cross-correlation  function: 

N 

A5^(r,x)  =  SJr.x)  -  5^(r,x)  s  ujT,t)  *  AM„(r,r)  .  (7) 

/i=0 

Using  first-order  perturbation  theory,  AS^(r,x)  may  be  related  linearly  taamplitude  and  dispersion 
perturbations  through  variations  to  ^e  amplitude  and  phase  kernels  of 


G„(K,r.t) 

=  G^(X,r,r)[l-»- Y„(X)] 

(8a) 

=  \(X,r,0  -  AGi„(k)t 

(8b) 

where  a  wavenumber-dependent  perturbation  to  the  relative  amplitude  and 

Aco  =  to  -  0)„  is  a  perturbation  to  the  dispersion  of  the  higher  mode.  In  terms  of  the 
differential  branch  cross-correlation  function, 


-oo 

A5^(r,t)  =  [C^(k,r,x)  yj})  +  B^(K,r,x)  A%(K)]  dk  (9) 

n  J  Q 

where  C^(X,r,x)  and  5^(X,r,x)  are  the  cross  correlations  of  the  mode-branch  synthetics  with  the 
Frdchet  mnels  for  the  amplitude  and  dispersion  perturbations.  This  linearization  is  valid,  provided 
that 
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In  general,  these  conditions  will  not  be  met,  and  the  inversion  must  be  iterated.  Equation  (9)  can  be 
rewritten  in  matrix  form. 


AS  =  Cy  +  BAo) . 


(11) 


By  applying  Rayleigh's  principle  and  first-order  perturbation  theory,  the  dispersion  can  be 
parameterized  in  terms  of  perturbations  to  the  reference  model  [Woodhouse  and  Dahlen,  1978]: 


Acd  =  HAm 


(12) 


where  H  is  the  Fnfchet  kernel  for  the  appropriate  model  parameter.  The  relationship  between  the 
amplitude  perturbation  and  the  model  perturbation  is  more  complicated,  however,  and  depends  on 
details  of  the  source  and  path.  Lerner-Lam  and  Jordan  [1983]  apply  an  orthogonalization 
procedure  to  reduce  the  sensitivity  of  the  branch  cross-correlation  functions  to  amplitude 
differences  betwwn  the  observed  and  synthetic  seismogram.  Their  approach  utilizes  a  projection 
operator  to  annihilate  the  matrix  C.  In  Ais  study,  we  consider  only  perturbations  to  the  dispersion 
and  disregard  perturbations  to  the  amplitude  kernel. 

With  this  formulation,  the  linearized  forward  problem  between  the  differential  branch 
cross-correlation  function,  AS,  and  perturbations  to  the  reference  earth  model.  Am,  reduces  to: 


AAm  =  AS  (13) 


where  A  =  BH.  So  far  we  have  not  considered  the  effect  of  noise  or  errors  in  the  data 
functionals.  In  general,  the  estimated  differential  branch  cross-correlation  function,  ASg,  will  be 
the  sum  of  the  true  branch  cross-correlation  function,  AS,  and  some  error,  n ,  such  that: 


ASg  =  AS  +  n  .  (14) 


We  assume  that  the  expected  value  of  n  is  zero. 


<n>  =  0  (1-^ 


and  therefore  its  covariance  matrix  is 


V„  =  <nnT>  .  (16) 


If  the  observational  errors  are  uncorrelated,  then  V-  will  be  diagonal.  This  is  almost  never  the 
case,  but  we  will  follow  in  the  footsteps  of  many  others  and  assume  it  to  be  true.  At  this  point, 
there  is  not  enough  known  about  the  error  processes  to  model  them  in  a  more  sophisticated  way. 

Rescaling  equation  13  to  account  for  the  observational  and  representational  errors: 


+  n  (17) 

where  Am  is  now  a  stochastic  process  whose  mean  satisfies  equation  13  and 

%  =  V^^'^A  (18a) 


'  ^  ^  1  I  V  ‘  ^  l  ^  1  "J  I  '  IV' ^  ' 


■^-TV'^ 


II 

c 

9  1 

(18b) 

A  -1/2 

n  =  V„  n  . 

(18c) 

With  this  scaling,  n  is  a  zero-mean  stochastic  process  whose  covariance  matrix  is  equal  to  the 
identity  operator  1.  In  practice,  we  assume  that  the  diagonal  elements  of  V  corresponding  to 
the lag- time  point,  r-, ,  of  the  i*  branch  cross-correlation  function  are  of  me  form  c^-  W(  r ,), 
where  W(r)  is  a  cosine4quare  taper  of  width  t.  Thus,  Oj-  is  the  standard  deviation  at  zero  lag. 

Lerner-Lam  and  Jordan  [1983]  solve  the  linearized  forward  problem  of  equation  17  using 
generalized  least-squares  and  imposing  the  conditions  of  regularization  and  smoothness; 


llA^g  -  ^Am  11  ^  -t-  a  11  Am  11^-1-  pllD^Amll^  =  min  (19) 

where  11  Ai^^  -  ^  Am  IP  is  the  least-squares  measure  of  data  misfit,  11  Am  IP  is  the  model 
perturbation  size,  is  the  second  difference  operator,  11  D^Am  I1  ^  is  the  model  perturbation 
smoothness,  and  a  and  p  are  the  scalar  trade-off  parameters  which  control  the  weighting  of  one 
term  relative  to  another.  The  regularization  has  the  effect  of  reducing  the  contribution  of  the 
poorly-determined  eigenvalues,  while  the  smoothness  constraint  controls  the  oscillations  of  the 
solution  about  the  starting  model.  Together,  they  stabilize  the  inversion.  The  solution  to  this 
minimization  problem  yields  an  estimate  of  the  model  perturbation.  Am 

Am  =  (20) 


where  is  an  operator  of  the  form; 


=  Ojk  al  -I-  PD2Td2)-1  P.  (21) 


Equations  20  and  21  constitute  the  generalized  least-squares  solution  to  the  forward  problem 
posed  by  equation  17.  This  formulation  of  the  waveform  inversion  technique  is  advantageous  in 
that  it  isolates  the  modes  in  the  time  domain,  obviating  the  need  for  complicated  processing  of  the 
data  in  the  frequency-wavenumber  domain.  In  particular,  this  approach  permits  the  contribution 
from  each  mode  to  be  assessed  and  weighted  individually.  Extensive  tests  with  synthetic  data 
[Lerner-Lam,  1982]  show  that  the  linearization  of  equations  8a  and  b  may  break  down  when  the 
starting  model  is  sufficiently  far  from  the  true  structure,  locking  the  inversion  into  a  spurious 
minimum.  However,  this  condition  may  be  identified  readily  by  direct  comparison  of  with  the 
synthetic  seismograms  with  the  data  and  corrected  by  initializing  the  inversion  with  an  appropriate 
perturbation  to  the  starting  model. 


Resolving  power  analysis 


Since  the  forward  problem  has  been  linearized  by  the  application  of  first-order  perturbation 
theory,  its  solutions  for  any  particular  data  set  may  be  evaluated  with  the  full  power  of  linear 
inverse  theory.  For  example,  we  can  compute  the  expected  value  of  the  estimated  model 
perturbation: 

(Am)  =  (  ^^A  §  +  n ) 


=  (^^  ^Am>  +  (^^n> 


= 


(22) 


is  called  the  model  resolution  operator  [Wiggins,  1972;  Menke,  1984]  and  Am  is  the 
exact  model  perturbation  which  satisfies  equation  13.  We  can  see  from  this  analysis  that  Rj^  acts 
as  a  filter  on  Am.  If  Rm  is  the  delta  function  which  satisfies  the  constraint  of  unimodularity,  then 
our  estimate  of  the  moael  perturbation  will  be  exact  In  general,  is  not  a  delta  function  and  our 
estimate  of  Am  is  a  weighted  average  of  the  true  model  perturbation.  A  solution  is  considered  to 
have  'good'  resolution  if  the  kernel  is  narrow  and  well  localized;  i.e.  if  it  is  peaked  and  has  small 
sidelobes. 

We  can  also  calculate  the  covariance  matrix  of  the  model  estimates,  Vj^: 


Vjjj  =  ((Am  -  (Am»  (Am  -  (Am))'^’) 


=  ^tT  , 


(23) 


Vj^  characterizes  the  errors  in  the  model  estimates  induced  by  errors  in  the  data.  The  diagonal 
elements  of  the  covariance  operator  represent  the  marginal  variance  of  the  model  values,  while  the 
off-diagonal  elements  measure  the  correlation  between  the  errors  at  one  depth  with  the  errors  at 
other  depths. 

The  properties  of  Rj^  and  depend  on  the  earth  model,  the  source  and  receiver  geometry, 
and  the  weighting  and  windowing  operators.  There  is  a  well-known  trade-off  between  resolution 
and  covariance  [Backus  and  Gilbert,  1970],  which  is  governed  by  the  parameters  a  and  /?in  this 
study.  Varying  a  and  ^  in  a  fixed  ratio  will  translate  the  solution  along  a  trade-off  curve  between 
the  ability  to  resolve  detail  and  the  reliability  of  the  estimate. 


EXPERIMENTS 

Although  Lerner-Lam  and  Jordan  [1983]  recognized  that  their  formulation  permitted  the 


resolving  power  analysis  by  standard  linear  methods,  they  did  not  use  it  to  assess  their  solutions. 
In  this  study,  we  incorporate  the  resolving  power  analysis  into  the  waveform  inversion  technique 
and  use  this  capability  to  examine  some  issues  in  experimental  design.  We  approach  this  problem 
by  designing  a  hypothetical  array  of  five  stations,  located  on  the  Baltic  Shield  and  distributed  over  a 
distance  of  1600  km  (Table  1).  This  receiver  distribution,  which  is  similar  to  the  configuration 
used  in  the  study  by  L^rner-Lam  and  Jordan  [1983],  is  modelled  after  the  NARS  array  in  order  to 
examine  the  contribution  of  such  dense,  linear  arrays  to  studies  of  upper-mantle  structure.  A 
cluster  of  five  earthquakes  which  occurred  in  the  two-year  period  1982-1984  is  selected  as  the 
source  array  (Table  2).  Located  south  of  Japan,  these  events  range  in  depth  from  1 19  to  552  km 
and  have  well-determined  centroids  and  moment  tensors  from  the  Harvard  solutions  [Dzicwonski  et 
ai,  1981, 1983, 1984].  Vertical-component  synthetic  seismograms  were  calculated  in  the  frequency 
band  from  5  to  35  mHz  using  Lemer-Lam's  asymptotic  travelling- wave  programs  with  Cara  et  al.'s 
[1980]  model  of  Eurasian  structure  and  the  attenuation  model  of  Masters  and  Gilbert  [1983]. 
Figure  2  is  a  map  of  the  source  and  receiver  distribution  considered  in  this  study  and  Figure  3  is  an 
example  of  the  synthetic  seismograms  generated  for  this  array  for  a  shallow  and  deep  event.  The 
base  station  of  this  array,  KONGO  1,  coincides  with  the  Global  Digital  Seismic  Network  station 
KONO.  Only  the  fundamental  and  the  first  four  higher-mode  branch  cross-correlation  functions 
were  included  in  the  resolving  power  calculations. 

With  this  configuration  of  sources  and  receivers,  we  explore  the  shear-wave  velocity 
resolving  power  of  the  higher-mode  waveform  inversion  technique.  In  particular,  we  consider  four 
end-member  cases;  a  single  source  recorded  at  one  and  five  receivers  and  five  sources  recorded  at 
one  and  five  receivers.  We  also  examine  the  influence  of  errors  in  the  data  by  experimenting  with 
the  relative  weighting  of  the  branch  cross-correlation  functions.  In  the  first  six  experiments,  we 
model  the  errors  as  stochastic  processes  by  weighting  the  branch  cross-correlation  functions 
equally.  This  'natural'  weighting  allows  the  mode  branch  with  the  highest  amplitude  (and 
consequently  the  highest  signal-to-noise  ratio)  to  control  the  inversion.  In  the  last  two  experiments, 
we  model  the  errors  as  processes  which  scale  with  signal  strength  by  weighting  the  branch 
cross-correlation  functions  accordingly.  All  the  figures  discussed  below  are  plotted  at  a  constant 
scale  to  facilitate  direct  comparison  of  the  results. 


Results 

In  our  first  two  experiments,  we  examine  the  resolving  power  of  a  single  source  recorded  at 
a  single  receiver  as  a  function  of  source  depth.  Figure  4  presents  the  resolution  and  covariance 
operators  for  the  source  at  1 19  km  depth.  These  operators  are  plotted  as  a  function  of  two  depth 
variables,  with  the  vertical  axis  identified  as  'target  depth.'  Each  row  of  may  be  characterized 
as  the  resolving  kernel  for  a  specific  target  depth.  For  example,  the  kernel  at  150  km  depth  ideally 
should  be  a  delta  function  centered  at  150  km.  Although  it  is  peaked,  this  particular  kernel  has  a 
width  of  approximately  200  km,  which  indicates  that  the  model-perturbation  estimate  at  150  depth 
is  a  weighted  average  of  the  true  model  perturbation  from  about  50  to  250  km.  The  absence  of 
sidelobes  signifies  that  the  averaging  is  localized.  The  amplitude  of  the  resolving  kernels  in  this 
experiment  falls  off  rapidly  with  depth  and  is  essentially  flat  by  300  km,  restricting  resolution  to  the 
top  of  the  upper  mantle.  The  covariance  operator  describes  how  errors  in  the  data  map  into  errors 
in  the  model  perturbation  estimate.  The  diagonal  elements  of  Vj^  represent  the  marginal  variance  of 
the  model  values,  while  the  off-diagonal  terms  measure  the  correlation  between  errors  at  one  depth 
with  errors  at  other  depths.  In  this  experiment,  the  variance  is  pronounced  at  shallow  depths,  with 
high  amplitudes  in  the  diagonal  and  off-diagonal  elements,  but  decreases  with  depth.  The  variance 
is  low  where  the  resolving  kernel  amplitude  is  small  and  higher  at  shallow  depths  where  the  kernel 
amplitude  is  large,  which  is  a  manifestation  of  the  trade-off  between  bias  and  covariance. 

Figure  5  depicts  the  same  experiment  for  a  single  deep  source  (/i=552  km).  In  contrast  to 
Figure  4,  the  resolving  kernels  are  peaked  down  to  target  depths  as  great  as  600  km.  These  kernels 
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are  generally  narrower  than  those  for  the  shallow  source,  although  the  width  increases  with  depth. 
Unlike  the  kernels  in  Figure  4,  however,  these  have  well-developed  sidebands,  which  indicate  that 
the  averaging  is  not  localized.  The  overall  increase  in  resolution  at  depth  is  mirrored  by  an  increase 
in  the  variance.  The  covariance  operator  for  this  experiment  has  high-amplitude  diagonal  terms  and 
considerable  structure  in  the  off-diagonal  elements. 

Comparison  of  Figures  4  and  5  illustrates  the  contribution  of  higher- mode  data  to  resolution 
of  upper-mantle  structure.  It  is  clear  from  Figure  1  that  the  seismogram  for  the  shallow  source  is 
dominated  by  the  fundamental  mode,  with  very  little  excitation  of  the  higher  modes.  With  the 
branch  cross-correlation  functions  all  weighted  equally,  the  fundamental  mode  controls  the 
inversion,  restricting  resolution  to  depths  shallower  than  200  km.  On  the  other  hand,  the 
seismogram  for  the  source  at  552  km  depth  contains  very  little  fundamental-mode  information;  the 
inversion  is  controlled  primarily  by  the  second,  third,  and  fourth  higher  modes.  Higher-mode 
dispersion  is  more  sensitive  to  structures  at  depth  than  fundamental-mode  dispersion,  and  this  is 
reflected  in  the  resolution  operator  of  Figure  5. 

In  the  next  two  experiments,  we  examine  the  effect  of  expanding  the  number  of  receivers 
from  one  to  five.  Figure  6  presents  the  results  from  the  experiment  with  the  deep  source  recorded 
at  the  array  of  receivers.  The  resolution  operator  for  this  experiment  is  smoother  and  the  peaks  are 
narrower  with  higher  amplitudes  than  the  kernels  in  Figure  5.  This  increase  in  resolution  is  not 
surprising,  as  the  array  spans  the  triplication  of  SSS,  a  phase  which  bottoms  between  400  and  800 
km.  Since  SSS  samples  the  lower  part  of  the  upper  mantle,  the  triplication  provides  valuable 
information  about  the  stucture  at  these  depths.  However,  the  greatest  gain  from  the  expansion  of 
the  receiver  array  is  observed  in  the  covariance  operator.  The  array  reduces  the  variance  of  the 
solution  by  nearly  a  factor  of  five,  since  we  assumed  that  the  errors  in  the  differential  branch 
cross-correlation  functions  are  uncorrelated.  Figure  7  displays  the  resolution  and  covariance 
operators  for  the  experiment  with  the  shallow  source  and  the  array  of  receivers.  Comparison  of 
Figures  7  and  4  shows  that  only  a  marginal  improvement  in  resolution  results  in  this  case,  although 
the  amplitude  of  the  covariance  operator  is  decreased  by  nearly  a  factor  •''f  five. 

Although  the  primary  benefit  from  an  array  of  receivers  appears  to  be  the  reduction  of  the 
covariance  operator,  we  can  take  advantage  of  the  relationship  between  resolution  and  covariance 
(which  is  governed  by  the  parameters  a  and  jS  in  equation  19)  to  exchange  some  of  the  decrease  in 
covariance  for  an  increase  in  resolution.  Figure  8  shows  the  results  from  increasing  ct  and 
decreasing  in  a  fixed  ratio  for  the  experiment  with  the  deep  source  recorded  at  the  array  of 
receivers.  Although  the  resolving  kernels  are  narrower  and  more  peaked  for  this  translation  along 
the  trade-off  curve,  the  cost  in  terms  of  the  reliability  of  the  solution  is  quite  high.  The  covariance 
operator  for  this  experiment  is  not  well  behaved,  with  errors  at  one  depth  scaling  with  errors  at 
nearly  all  depths.  Figure  9  presents  the  results  for  a  translation  along  the  trade-off  curve  in  the 
other  direction.  In  this  case,  the  resolution  operator  is  highly  degraded,  while  the  covariance 
operator  is  nearly  flat.  We  can  also  experiment  with  the  trade-off  parameters  for  the  shallow  source 
recorded  at  the  array  of  receivers,  resolution  is  limited  more  by  the  depth  of  the  source  (given  the 
'natural'  weighting  of  the  branch  cross-correlation  functions)  than  by  the  number  of  receivers  in 
this  experiment. 

Having  examined  the  dependence  of  the  resolution  and  covariance  operators  on  an  array  of 
receivers,  we  model  the  effect  of  an  array  of  sources.  Figure  10  shows  the  results  from  an 
experiment  with  five  sources  distributed  in  depth  from  119  to  552  km  and  recorded  at  a  single 
station.  From  a  comparison  of  Figures  4  and  5  with  Figure  10,  it  is  clear  that  a  distribution  of 
sources  can  be  particularly  effective  in  enhancing  resolution.  The  resolution  kernels  of  Figure  10 
have  high-amplitude  and  narrow  peaks  and  show  some  resolution  as  deep  as  600-700  km  depth. 
The  sideband  structure  is  damped,  although  the  averaging  is  not  completely  localized.  In  terms  of 
the  covariance  operator,  an  array  of  sources  reduces  the  variance  considerably,  but  not  as 
effectively  as  an  array  of  receivers.  Figure  1 1  illustrates  the  resolution  and  covariance  operators  for 
the  optimal  experiment  of  five  sources  recorded  at  five  receivers.  As  before,  the  addition  of  five 
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receivers  narrows  the  peaks  of  the  resolution  kernels,  but  the  greatest  gain  is  in  the  reduction  of  the 
variance. 

In  these  six  experiments,  the  branch  cross-correlation  functions  are  effectively  weighted  by 
their  relative  amplitudes  in  the  inversion,  since  the  error  at  zero  lag,  ,  is  assumed  to  be 
independent  of  the  branch  number.  This  'natural'  weighting  entails  certain  implicit  assumptions 
about  the  error  process.  If  one  assumes  that  errors  in  the  data  are  introduced  strictly  by  the 
presence  of  ambient  seismic  noise  (wind,  microseisms,  cultural  noise,  etc.),  then  this  natural 
weighting  by  the  relative  excitation  of  the  modes  is  most  appropriate.  In  that  case,  the  highest 
amplitudes  carry  the  greatest  weight  since  they  have  the  best  signal-to-noise  ratio.  However,  many 
errors  may  be  a  function  of  signal  strength  (such  as  multipathing,  source  structure,  etc.).  If  these 
errors  dominate  the  random  processes,  then  the  natural  weighting  is  not  appropriate. 

We  explore  the  effect  of  characterizing  the  error  processes  by  experimenting  with  the  mode 
branch  weights  for  the  shallow  source  recorded  at  one  and  five  receivers.  By  increasing  the  weight 
of  the  higher-mode  branches  relative  to  the  fundamental,  we  model  the  case  where  the  errors  scale 
with  signal  strength.  Figure  12  illustrates  the  resolution  and  covariance  operators  for  this 
experiment.  Comparing  Figures  4  and  12,  it  is  clear  that  this  weighting  of  the  higher-mode  signals 
has  improved  the  resolution  at  depth,  although  these  kernels  are  neither  narrow  nor  localized.  The 
increase  in  resolution  is  reflected  by  an  increase  in  the  variance.  As  before,  increasing  the  number 
of  receivers  (Figure  13)  marginally  improves  the  resolution,  but  dramatically  decreases  the  variance 
of  the  solution.  The  results  from  these  two  experiments  highlight  the  importance  of  the 
assumptions  made  about  the  error  process.  If  errors  scale  with  signal  strength  (which  we  think  is 
more  realistic),  then  the  branch  cross-correlation  functions  need  to  be  weighted  accordingly. 
However,  the  weighting  of  the  branch  cross-correlation  functions  plays  an  importani  role  in  the 
resolving  power  of  a  particular  experiment,  as  it  is  possible  to  obtain  some  resolution  at  depth  even 
from  a  shallow  source  when  the  contributions  from  the  higher-mode  branches  are  weighted 
appropriately. 


Conclusions 

In  this  study,  we  augment  the  waveform  inversion  technique  of  Lerner-Lam  and  Jordan 
[1983]  with  a  resolving  power  analysis  and  use  this  capability  to  examine  some  issues  in 
experimental  design.  These  experiments  have  demonstrated  several  important  results.  First, 
fundamental-mode  data  alone  cannot  resolve  structure  at  depths  greater  than  200  km;  a  fact  which 
has  been  appreciated  by  seismologists  for  many  years.  Higher-mode  information,  which  samples 
greater  depths,  is  essential  to  obtain  resolution  of  the  upper  mantle.  Second,  a  single  station  is 
capable  of  nearly  the  same  resolving  power  as  a  dense  array  of  receivers  with  this  method,  since  the 
higher-mode  waveform  inversion  technique  does  not  rely  on  spatial  transforms.  The  advantage  of 
such  an  array  comes  primarily  through  the  reduction  of  the  variance  of  the  solution,  although  the 
decrease  in  covariance  may  be  exchanged  for  an  increase  in  resolution  through  the  trade-off 
parameters.  Third,  source  arrays  are  particularly  effective  in  enhancing  resolution  at  depth, 
provided  that  the  source  depths  are  well  distributed  and  the  source  centroids  and  moment  tensors 
are  well  determined.  Finally,  the  characterization  of  the  error  processes  plays  an  important  role  in 
the  resolving  power  analysis  of  a  particular  experiment,  since  the  relative  weighting  of  the  branch 
cross-correlation  functions  will  influence  the  resolution.  Since  the  waveform  inversion  technique 
permits  the  contribution  of  each  branch  cross-coiTelation  function  to  be  assessed  and  weighted 
separately,  resolution  of  structure  at  depth  is  possible  even  with  shallow  sources. 
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Station  I.D, 

Latitude  (°) 

Longitude  (°) 
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K0N017 

67.989 

41.243 

Date 

Latitude(°) 

Table  2 

The  Source  Arrav 

Longitude  (°) 

Depth  (km) 

Magnitude 

07/04/82 

27.92 

136.48 

551.8 

6.3 

07/05/82 

30.77 

130.47 

119.0 

5.7 

09/06/82 

29.18 

140.65 

155.6 

6.6 

01/01/84 

33.38 

136.81 

383.6 

6.5 

02/13/84 

25.26 

122.29 

268.2 

5.5 

23 


TIME  (min)  TIME  (min) 

Figure  3:  This  figure  illustrates  the  synthetic  seismograms  calculated  for  the  receiver  arri 
the  events  07/04/82  and  07/05/82.  The  stations  span  the  triplication  of  SSS,  which  prov 
critical  information  about  the  lower  part  of  the  upper  mantle. 
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RESOLUTION  COVARIANCE 

SOURCE  (H«552  KM).  5  RECEIVERS  1  SOURCE  (H«552  KM).  5  RECEIVERS 


Figure  6:  This  flgurc  displays  the  resolution  and  covariance  operators  for  a  deep  source  recorded  at 
the  array  of  receivers.  The  resolving  kernels  have  been  smoothed  and  their  peaks  narrowed  in 
this  experiment,  while  the  covariance  has  come  down  by  a  nearly  a  factor  of  five. 


Figure  7:  The  resolution  and  covariance  operators  for  the  shallow  event  recorded  at  five  receivers 
are  plotted  in  this  figure.  For  constant  values  of  a  and  P,  the  addition  of  a  receiver  array  results 
in  marginal  improvement  in  the  resolution,  although  the  covariance  is  cut  down  by  nearly  a  factor 
of  five.  In  this  experiment,  the  resolving  power  is  constrained  more  by  the  depth  of  the  source 
than  by  the  number  of  receivers. 


Figure  8:  Since  the  receiver  array  has  reduced  the  covariance  so  dramatically,  it  is  possible  to 
experiment  with  the  trade-off  curve  and  exchange  the  reduction  in  covariance  fcM-  an  increase  in 
resolution.  Figure  8  demonstrates  the  results  from  an  inversion  where  a  andp  were  varied  in  a 
fixed  ratio.  Comparing  Figures  6  and  8,  it  is  clear  that  the  resolution  can  be  improved,  although 
the  cost  in  terms  of  the  increase  in  the  covariance  is  quite  high. 
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Figure  11:  This  figure  presents  the  resolution  and  covariance  operators  for  the  optimal  experiment 
of  five  sources  recorded  at  five  receivers.  As  before,  the  addition  of  the  receiver  array  has 
smoothed  the  resolving  kernels  and  narrowed  the  pe^s,  while  damping  the  covariance.  Similar 
to  experiment  of  Figure  8,  this  reduction  in  the  covariance  may  be  trad^  for  an  increase  in 
resolution. 


Figure  13:  This  figure  expands  the  noise  model  experiment  to  the  shallow  source  recorded  at  five 
receivers.  As  before,  the  addition  of  the  receiver  array  gives  a  reduction  in  covariance  which 
may  be  traded-off  for  an  improvement  in  resolution. 


